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Abstract 


A parallel, finite-volume algorithm has been 
developed for large-eddy simulation (LES) of 
compressible turbulent flows. This algorithm includes 
piecewise linear least-square reconstruction, trilinear 
finite-element interpolation, Roe flux-difference 
splitting, and second-order MacCormack time 
marching. Parallel implementation is done using the 
message-passing programming model. In this paper, the 
numerical algorithm is described. To validate the 
numerical method for turbulence simulation, LES of 
fully developed turbulent flow in a square duct is 
performed for a Reynolds number of 320 based on the 
average friction velocity and the hydraulic diameter of 
the duct. Direct numerical simulation (DNS) results are 
available for this test case, and the accuracy of this 
algorithm for turbulence simulations can be ascertained 
by comparing the LES solutions with the DNS results. 
The effects of grid resolution, upwind numerical 
dissipation, and subgrid-scale dissipation on the 
accuracy of the LES are examined. Comparison with 
DNS results shows that the standard Roe flux-difference 
splitting dissipation adversely affects the accuracy of the 
turbulence simulation. For accurate turbulence 
simulations, only 3-5 percent of the standard Roe flux- 
difference splitting dissipation is needed. 

Nomenclature 


A area 

\A\ Roe flux-difference splitting matrix 

A s area of duct side walls 

c local speed of sound 

C SGS model constant 
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CFD 

computational fluid dynamics 

CFL 

Courant number 

c s 

Smagorinsky constant 

c v 

specific heat at constant volume 

d 

normal distance from a solid wall 

ct 

normal distance from a solid wall in wall 

J+ P U X d 
units, d = 

D 

entire flow domain 


hydraulic diameter 

DNS 

direct numerical simulation 

(6 

elemental surface area on the boundary of a 
control volume 

dv 

elemental volume of a control volume 


total energy/unit volume 

f 

normal component of the inviscid flux 
vector 

t 

flux vector of the Navier-Stokes equations 

FDS 

flux-difference splitting 


inviscid flux vector 


viscous flux vector 

G 

spatial filter used in the LES equations 

I 

total number of grid points in the 
streamwise direction 

J 

Jacobian determinant 

J 

total number of grid points in the wall- 
normal direction 

k 

conduction heat-transfer coefficient 

K 

total number of grid points in the spanwise 
direction 

LES 

large-eddy simulation 

MPI 

Message-Passing Interface 
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a 

normal unit vector 

Subscripts 


p 

static pressure 

a 

node index 

P s 

mean pressure gradient 

L 

flow conditions to the left of a cell face 

PVM 

Parallel Virtual Machine 

rms 

root mean square 

2 

q 

trace of the SGS Reynolds stress tensor 

R 

flow conditions to the right of a cell face 


SGS term in the LES energy equation 

Superscript 


R 

specific gas constant 

n 

time level 

S 

total area on the boundary of a control 
volume 


cell-averaged quantities in the Navier- 
Stokes equations, or filtered or large- 

SGS 

subgrid scale 


scale quantities in the LES equations 


velocity gradient tensor 


Favre- filtered (density-weighted) variables 


t time vector quantity 


T 

temperature 

u 

x-component velocity 

u ave 

mean streamwise velocity 

~«v 

2 

mean Reynolds stress 



T 

W T 

friction velocity, w T = 

u 

state vector of the Navier-Stokes equations 

V 

y-component velocity 

V 

total volume of a control volume 

w 

z-component velocity 

x, y, z 

coordinates of the physical space 


Kronecker delta 

A 

width of filter used in LES equations 

A ', 

sampling time 

e l 

scaling factor for Roe flux-difference 
splitting 

P 

molecular viscosity coefficient 

5.T1.S 

coordinates of the computational space 

P 

static density 


SGS term in the LES momentum equations 
(the SGS Reynolds stress tensor) 

X kl 

viscous stress tensor 

V 

wall shear stress 


Introduction 

Turbulence dominates the internal flows in aircraft jet 
engine components such as inlets, ducts, and nozzles 
and has been found to significantly influence engine 
noise and performance. Analytical tools are therefore 
needed to provide accurate predictions of these 
turbulent f ows and allow engineers to explore the 
underlying flow physics, which would allow better 
aeropropuhion flow components to be designed and 
used in the aerospace industry. 

Direct numerical simulation (DNS) of the turbulent 
flow inside of complete jet engines is presently not 
possible because of the tremendous computational 
resources required; however, technologies that 
potentially :ould make such a feat possible in the future 
are availab e today. These technologies include large- 
eddy simuhtion (LES) of turbulent flows, unstructured 
computatio lal fluid dynamics (CFD) algorithms, and 
parallel computer systems. Large-eddy simulation has 
been shovm to provide accurate turbulent flow 
simulation it a fraction of the cost of direct simulation. 
With unstiuctured CFD algorithms, complex three- 
dimensional aerodynamics shapes, including complete 
aircraft gee metrics, have been modeled using a single 
grid. Largo-eddy simulation and unstructured CFD 
algorithms require large computing resources that 
potentially can be provided by the emerging parallel 
computer systems. By linking together hundreds or even 
thousands of individual processor nodes, the parallel 
computer systems can deliver significant advances in 
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computational resources in terms of memory, storage, 
and computing speed. 

The above three technologies have been the subjects 
of ongoing intensive research, and a large body of 
knowledge has been separately accumulated on each of 
these subjects. The objective of this research is to 
develop a turbulence simulation tool using a 
combination of all of these technologies. The accuracy 
and efficiency of such a tool for turbulence simulations 
are then examined in detail from the LES of fully 
developed turbulent flow in a square duct. 

Use of trade names or names of manufacturers in this 
document does not constitute an official endorsement of 
such products or manufacturers, either expressed or 
implied, by the National Aeronautics and Space 
Administration. 

Numerical Algorithm 

Development of the numerical algorithm has 
previously been described in detail. 1 This algorithm has 
previously been validated for time-accurate inviscid 
Euler simulations 2 and three-dimensional viscous 
Navier-Stokes simulations with good results. To 
describe the numerical algorithm, the Navier-Stokes 
equations are used in this section. These equations can 
be written in vector form as 

^ + ^•£ = 0 (i) 

dt 

where U is the state vector and P is the flux vector of 
the Navier-Stokes equations. 

The above equation is discretized using the finite- 
volume approach. In this approach, equation (1) is 
integrated over a finite volume. Assuming the grid does 
not change with time and using the Gauss divergence 
theorem, the resulting equation is 

jp = ® 

where 

u = pju dv (3) 

v V 


and 

Pd* = f tids (4) 

To numerically solve equation (2), the major steps of 
the solution procedure are reconstruction, flux 
computation, and evolution. This standard, finite- 
volume solution procedure has been used in previous 
works and has been described in detail by Barth 3 4 The 
steps are given below. 

Step One: Reconstruction 

For the first step, reconstruction, a cell-centered 
scheme is used. The piecewise linear, least-square 
reconstruction procedure used here is similar to those 
used by Barth 4 and Coirier. 5 Each of the five primitive 
variables p , w, v, w, and p is assumed to linearly vary 
within a finite volume as: 

U(x, y, z) = 0 + U (x-x) + U (y-y) 

y (5) 

+ U z (z-z) 

where U can be any of the above variables. The bars in 
equation (5) denote cell-averaged values as defined in 
equation (3). When used for high-speed compressible 
flow simulations, a gradient limiter is normally used in 
equation (5) to ensure that the reconstruction 
polynomial does not produce new extrema near a flow 
discontinuity such as a shock wave. In this paper, the 
gradient limiter is not used because the test case is low- 
Mach number turbulent flow in a square duct. 

Following Coirier, 5 the gradients U x , U , and U z in 
the target cell are computed using a least-square 
procedure that minimizes the sum of the squares of the 
differences between the values computed using the 
reconstruction polynomial from the target cell and the 
cell averages of the support set. For a three-dimensional 
hexahedron cell, the support set is the six neighboring 
cells that share their faces with the target cell. 
Algebraically, the minimization statement above can be 
expressed as: 


b j b 0 


U x 


C 1 

b | a 1 b 3 


U y 

- 

c 2 

bry 1?^ 


U z 


-°3 


3 

American Institute of Aeronautics and Astronautics 



6 

«i = X (\-*o > 2 

i = 1 
6 

a 2 = Z(y,-y 0 ) 2 

i = 1 
6 

«3 = K*,-*0> 2 

I = 1 
6 

b \ = X (*,■- W.-yo) 

i = 1 
6 

= X (^ -*o)(z. - 2 o) ( 7) 

I = 1 
6 

b i= X (y,-yo)( 2 .- 2 o) 

i = I 
6 

c, = X (O, -O o )(x ( .-x 0 ) 

I = 1 
6 

c 2 = I (O, ~O o )(y,-y 0 ) 

i = i 
6 

*3 = X (U / ~U 0 )(Z | -Z 0 ) 

i = 1 

where i = 1-6 denotes the neighboring cells, and i = 0 
denotes the target cell. 

Step Two: Flux Computation 

With the piecewise linear reconstruction, the 
unknown variables are continuous and assumed to 
linearly vary within a finite volume. However, no 
guarantee exists that the variables will be continuous 
across adjacent finite volumes because a different 
polynomial is used in each finite volume. As a result, a 
flux formula is needed to compute a single flux at a 
finite- volume boundary using fluxes from the adjacent 
volumes. In the numerical solution of the Navier-Stokes 
equations, splitting the total flux vector into the inviscid 
flux vector and viscous flux vector is convenient: 

^ ( 8 ) 

For the viscous flux, a simple arithmetic average is 
used. The normal component of the inviscid flux 
vector, rt , is approximated using the Roe flux- 
difference splitting (FDS) method without the entropy 
correction. The entropy correction is normally used to 
remove the nonphysical expansion shock at the sonic 


transition p^int and is not needed here because of the 
low-Mach lumber test cases. 

Define tl e normal component of the inviscid flux 
vector as 


f = r, • a ( 9 ) 

Then f can be computed using the Roe FDS method 
as 

f = |( f L + f R)-^ A l( U R- U L) (1°) 

where the subscripts L and R denote the flow conditions 
to the left and right sides of the cell face. 

Figure 1 shows the definitions used for the left and 
right states Consider a cell, A, and its neighbor, B, 
sharing a common face, 1-2. When the flux across face 
1-2 is computed for cell A, the left state (L) of the face 
1-2 is on the side of cell A, and the right state (R) is on 
the side of cell B. This definition of the left and right 
states of a face is used because the face normal unit 
vector fi, which also serves as the locally one- 
dimensional coordinate system for the wave 
propagation across face 1-2, points from cell A to cell 
B. The L and R states are reversed when the flux is 
computed for cell B. 



Figure 1. Definition of the left and right states of a face. 


Step Three: Evolution 

The twc-stage, second-order, MacCormack time- 
marching algorithm is used to advance the solution in 
time. This explicit predictor- corrector time-marching 
method is accurate, efficient, and simple to implement 
on parallel romputer systems. 

Equation (2) can be rewritten as 

J, B ‘ R (ll > 
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When applied to equation (11), the MacCormack time- 
marching method gives 

0" + ' = 0" + A/R n (12) 

U" +I = ^C" +1 + U n + AfR n + ‘j (13) 

In addition to the major solution steps outlined above, 
describing how the volume and surface integrals in 
equations (2) and (3) are evaluated is important. 
Although the flow variables are approximated by 
discontinuous piecewise linear polynomials, the spatial 
coordinates x, y, and z of a finite volume are 
approximated by a continuous trilinear hexahedral 
element. 5 6 This approach is the same that finite-element 
methods use to approximate the spatial coordinates. All 
of the integrations are then numerically evaluated using 
the one-point Gauss quadrature formula. 

Each cell in the physical £(x, y, z) space is mapped 
to a trilinear hexahedral element in the T|, £) space 
as shown in figure 2. The nodes, indexed 1 to 8, have the 
nodal coordinates in the ^ space shown in table 1. 



Figure 2. Local mapping between the physical finite- 
volume and the trilinear hexahedral element. 


Table 1. Nodal coordinates in the q space. 


Node index (a) 

*>a 



1 

-1 

-1 

-1 

2 

1 

-1 

-1 

3 

1 

1 

-1 

4 

-1 

1 

-1 

5 

-1 

-1 

1 

6 

1 

-1 

1 

7 

1 

1 

1 

8 

-1 

1 

1 


This type of mapping is different from the mapping 
used in generalized curvilinear finite-difference 
methods. In the finite-difference methods, the mapping 
applies to the entire computational block. In this 
algorithm, the local mapping applies to the local cell 
only. Each cell has its own mapping function, which is 
similar to a finite-difference multiblock method in 
which each cell is its own block. This ability gives 
considerable flexibility in the grid topology that can be 
used and is one of the strengths of this unstructured 
finite-volume algorithm. 

The mapping from the | space to the & space is 
given by: 

4) = X N a& x a 

a = 1 
S 

4) = X N a&y a < 14 > 

a = 1 
8 

4) = X N a&K 

a = 1 

N a (b = i(l+^)(l+il a Tl)(l+; fl O (15) 

where the subscript a denotes the node index, ranging 
from 1 to 8. In the physical (x, y, z) coordinate system, 
node a has the coordinate (x , y , z ). In the 
computational (£, T|, £) space, node a has the 
coordinate (£, a ,r\ a , C, a ) . The coordinates (x u> y fl , z ) 
vary from cell to cell, depending on the physical grid. 
The coordinates r\ a , £ a ) are the same for every cell 
and are shown in table 1 . 

To evaluate the volume integral, the following 
relation 6 is used: 

J /(x, y, z )dv 

v 

= J J | /(x(£, Tl, 0, y(^, Tl, Q, z(£, T|, 0) (16) 

- 1 - 1-1 

j&T 1. Wnd? 

where j is the Jacobian determinant, defined as 



det 


\ x n x ; 

H U y; 


L z ^ z n z d 


(17) 
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The partial derivatives x^ , y^ , , and so forth can be 

obtained by differentiating equation (14). Evaluating the 
determinant in equation (17) results in the following: 


J = x ^ z ri" X 5 y ? Z h" X ;V5 + VcH 
+ x *Vt“ V5 Z <; 


(18) 


developed below. Approximations for the y and z 
coordinates are made exactly the same way. 


x 



(23) 


To evaluate the integral in equation (16), the one* 
point Gauss quadrature formula is used. In one 
dimension, the Gauss quadrature formulas are optimal, 
which means that accuracy of order (2 n) is achieved 
using (n) integration points. Gaussian rules for integrals 
in several dimensions are constructed by employing the 
one-dimensional Gaussian rules on each coordinate 
separately. In three dimensions, the one-point Gaussian 
rule is given as 

l 1 1 

J J j/s.% 

-1-1-1 ^ ' 

= 8 /ft = 0, H = 0, C = 0) 


i l l 

111 X&Tl.Oj&Tl 


\ | Jj&Ti ,Qdt,dr\dt 

-1-1-1 

_ 8x(^ ^ Q,rj = 0 , g = 0)j(£ = 0,n = 0, C = 0) (7 ~ 

8j(^ = 0,ti = 0,C = 0) 
or 

x = x(£ = 0, ti = 0, £ = 0) 

y = y(^ = 0,Ti = 0,C = 0) ( 26 ) 

z = z(£ = 0, r| = 0, £ = 0) 


Using the tools described above, the volume of the cell 
is computed as 

V = jldv (20) 

V 

l 1 1 

v = j J (2D 

-1-1-1 

Using equation (19), 

V = 8 j(£ = 0, Tl = 0, £ = 0) (22) 

so that the cell volume is approximately eight times the 
Jacobian determinant evaluated at the center of the cell 
£ = 0, r| = 0, £ = 0. Note that equation (22) contains 
two approximations: the physical coordinates (x, y, z) in 
the cell are approximated by equation (14), and the one- 
point Gaussian rule given by equation (19) is used for 
the numerical integration. Better approximation of the 
cell volume can be obtained using a higher-order 
approximation for the physical coordinates and a 
Gaussian rule with more points. 

Computing the centroid of the cell also requires the 
evaluation of the volume integral. The coordinates of a 
cell centroid are given by (x, y, z). The numerical 
approximation for the x coordinate of the centroid is 


From equation (14), 

. 8 

x(E = 0, Ti = 0, ; = 0) = ^ I x a 

a = 1 

l 8 

y(E = 0,Tl = <U = 0) = i X y a < 27 > 

a = 1 
8 

z(£ = 0,11 = 0, 5 = 0) = J X z fl 

a = 1 

With the approximations in this algorithm, the 
coordinates of the cell centroid are simply the averages 
of the coordinates of the eight nodes defining the three- 
dimensional hexahedron cell. 

The tasl of evaluating the surface integrals is 
described n ?xt. The surface integrals on the right side of 

equation (7) are of the form jt dS. To develop the 

procedure, >ne face of the hexahedron shown in figure 2 
is considered. The results for the other faces can be 
obtained us ng the same procedure. 

Consider face 6-2-3-7 of the hexahedron in figure 2. 
Figure 3 shows this face redrawn for convenience. 
Because ^ = 1.0 for the nodes 6, 2, 3, 7 as well as any 
other point on this face, equation (14) reduces to the 
following: 
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X 


(32) 


= 4 £ o+viXi+^K 

<1 = 6,2, 3,7 

y = ! S (i + vDO + C a C)y a (28) 

a =6,2, 3,7 

z = z X (i +voo + ; a 0z a 


di = di^ xdi x 


a = 6, 2, 3, 7 

so that x, y, and z are functions of T| and £ only. 



Figure 3. Coordinate systems for surface integral 
evaluation. 


Following the development outlined in Greenberg, 7 
the tangent vectors di j and di 2 are tangents on the 
plane 6-2-3-7. di x is defined to be along the 
r| = constant curve on the face, and di^ is along the 
£ = constant curve. The vector di x may be expressed as 
follows: 


di j = dx\ + dy] + dzk 
dx = x^d^ + x^dr\ + x^dt^ 
dy = y^d% + y^dr\ + 
c/z = z^d% + z^z/ri + z^d£ 


(29) 


where = 0 because the entire plane 6-2-3-7 is an 
i; -constant plane, and dv\ = 0 because the vector di x 
is defined to be along the T| = constant curve. So 


= (x^ + yj+z^lc)^ 

and similarly 

d$ 2 = (x^i + y^j + z n ic)dTi 


(30) 


(31) 


The elemental area vector di , denoted by the shaded 
parallelogram, can be computed as 


or 


d * = (yn z ?~Vd' + ( V;~Vs )J 

+ (x„y 5 -y^x 5 )i] dT]dC, 


(33) 


Note that the order of the cross product in equation (32) 
is chosen so that the elemental area vector di is positive 
pointing out of the cell and negative pointing into the 
cell. With the vector t defined as 


t = F x i + F y j + F z lc 


(34) 


the dot product is 

tdi = [(y T1 z^-z tl y ; )F x + (z T1 x ? -x T1 z c )F y 

+ (x„y 5 -y^x ? )F z ]rfrirfC 


(35) 


Finally, using the one-point Gaussian rule, the surface 
integral can be evaluated as 


j $d$ = 4[(y^z ? -z n y ; )F x 

(6237) 

+ <VrVc )F y 
+ ( V?" Vd F d 


(36) 

5= i.tjso.c-0 


The surface integrals for the other faces can be 
approximated in an analogous fashion. The results are 
given below. 

J P-dS = 4[(z T1 y i .-y 11 z !; )F x 

(1584) 

+ ( x n z C - VC )F y (37) 

+ (y^- x ^)F z ]| 4= _ 1>T1=0 , ;=0 


| # di = 4[(z^y^-y^z^)F x 
(8734) 

+ (x^-z^XjOFy 

+ (y^-x^)FJ 
j Pd$ = 4[(y 5 z ? -z 5 y ; )F x 

(1265) 

+ (z^-x^)F y 


(38) 


4 = o,ti = 1,5 = o 


(39) 
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j f -«ft = 4[(y 5 z 11 -z l jy 11 )F x 

(5678) 

+ ( z ^ x r| _ x ^ z r|)Fy 

+ U^-y^FJ 


(40) 


5 = o, ti = o, ; = 1 


3 £7 3 - - 3 3-- 

T, *d7 l E ' u ' + 17 l ’’“i--S7 l “^‘ 


3 ( i 3 f \ . do*, 

- 5T,r37,T “*-3^ 


*/ , dg, 

+ ° v 3x^ = ° 


(45) 


J tdl = 4[(z^y T1 -y^z T1 )F x 

(2143) 

+ ( x ^-HN )F y (41) 

*<V„-, 5 y„>F,l| 1 ,_o, 

Symbolically, the surface integrals of the opposite 
faces are the negative of each other (for example, face 
£ = 1 as given by equation (36) and face E, = -1 as 
given by equation (37)). However, for actual numerical 
values, each of the integrals will need to be separately 
evaluated because the integrands depend on the 
coordinates of the nodal points (x fl , y fl , z Q ) and the 
Gaussian point coordinates. 


The bar in the LES equations (43) to (45) denotes a 
filtered or large-scale flow quantity, defined as 

/ = jc(* -*')/(* vr 

D 

where G is a spatial filter and the integral is over the 
flow domain, D. The tilde in the LES equations denotes 
a Favre-filtered (density-weighted) variable, defined as 



The filtered ideal gas equation of state is given by 


To calculate the Jacobian determinant of the 
coordinate transformation, equations (14) and (15) are 
used. For example, the derivative x^ can be evaluated as 
follows: 


X(t)= I Njbx, 

a = 1 
8 

X $ = I N a,^a 

a = 1 


(42) 


Governing Equations for Large-Eddv Simulation 


In LES, the large scale of turbulence is computed 
directly in the numerical simulation, and the effects of 
the small scale stresses are modeled using a subgrid- 
scale (SGS) model. The governing equations for LES of 
turbulent flows can be obtained from filtering (local 
volume-averaging) the compressible Navier-Stokes 
equations. From Moin et al., 8 the LES equations for 
compressible flows (using tensor notation) are given by 
the following: 


a. a — 

5?+j^P“» = ° 


(43) 


3 3 3 ^ki 

37 pM * + 3 ^*“' + 3^-4 + 4 ° W) 


p = pR T 

Moin et al. 8 used an internal energy equation in their 
derivation of the LES equations. The LES total energy 
equation, equation (45), is obtained from adding the dot 
product of he LES momentum equation, equation (44), 
and the filtered velocity field u k to the LES internal 
energy equation in Moin et al. 8 

The LES equations given by equations (43) to (45) are 
essentially the Navier-Stokes equations written for the 
filtered variables plus the additional subgrid terms in the 
momentum and total energy equations. Thus, the 
numerical 2 lgorithm developed in the last section can be 
used to sol /e the LES equations. The treatment of the 
subgrid ten is are discussed in the next section. 

Subgrid-Scale Models 

Detailed studies have previously been performed to 
assess the relative importance of the subgrid terms in the 
filtered total energy equation for compressible turbulent 
shear flow; at different Mach numbers. 9, 10 These 
studies led to the conclusion that the energy subgrid 
terms may be neglected if the Mach number of the 
simulation s low. Because of the low Mach number of 
the turbule it square duct test case, this assumption is 
used. 
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The subgrid term in the momentum equations, 
equation (44), is 

°ki = p 

To close the system of LES equations, this term needs to 
be modeled. In Moin et al., 8 this term is approximated 
as follows: 

a kl = - 2CpA 2 * |s|js t/ - \s mm & k ^ + ] -q 2 8 kl (46) 

where 



is the trace of the SGS Reynolds stress tensor. The 
filtered velocity gradient tensor is 

1 (du k duA 
S kl = ^ ^ + 

2 t 3x / 3x J 

and 

l 

\s\ = (2 S k ,S kl ) 2 

In equation (46), C is a constant to be determined 
according to the particular SGS model used. For LES of 
turbulent channel and duct flows using the Smagorinsky 
SGS model, 11 a value of C = 0.01 is commonly used 
with good results. Note that the constant C in 
equation (39) is the square of the Smagorinsky constant 

c s =0.1. 

Unlike LES of isotropic turbulence, C is not constant 
in wall-bounded flows and varies according to distance 
from the wall. The dynamic SGS model developed by 
Germano et al. 12 would correctly determine the value 
for C using a dynamic procedure; however, this model is 
computationally expensive because of the extra filtering 
operations that must be done. Also, a question currently 
exists on the mathematical well-posedness of this 
model. Finally, the dynamic model has been known to 
compensate for the effects of numerical dissipation by 
automatically varying the magnitude of the constant C. 

One of the main objectives of this research is to 
quantify the effects of the upwind numerical dissipation 
on the accuracy of the turbulence simulations. Also, the 
simpler Smagorinsky model has been found to work as 
well as the dynamic SGS model for this simple test 
case. 14 As the result, the Smagorinsky SGS model is 
used in this study with the constant C given by the 
following: 


C = 0.0l(j.0-exp(-(j^ J) (47) 

where d + is the normal distance from the wall in wall 
units, defined as 



and the friction velocity is defined as 



Because the turbulent flow in the comer of the square 
duct encounters walls in two different directions, d is 
taken to be 

d = 2yz _ (48) 

2 o 2 
y + z + (y‘ + z“) 

Equation (48) is frequently used in the turbulence 
modeling of flows in the vicinity of a wall comer. The 
variables y and z are the normal distances to the nearest 
walls in the y and z directions. Note that d tends to y as 
(y/z) tends to 0, and d tends to z as (z/y) tends to 0, 
which are the intended results. 

In LES, the width of the filter used in the process of 
volume-averaging the Navier-Stokes equations, A , is 
typically chosen to be the grid spacing size. This study 
defined the grid spacing size to be l/V , where V is the 
cell volume. 

The term q 2 in equation (46) is the isotropic part of 
the SGS Reynolds stress tensor. Like the rest of this 
tensor, the term cannot be calculated directly in an 
LES and has to be modeled. A number of different 
models for q 2 has been proposed. 15 " 17 However, 
results from recent studies indicated that this term is 
not important for accurate LES of low-Mach number, 
low-Reynolds number compressible turbulent flows. 

An evidence in support of the above conclusion was 
presented by Squires, 18 who compared two different 
models of q 2 in addition to setting q 2 = 0. Squires found 
essentially no difference in the results of LES of 
compressible isotropic turbulence at a low Mach 
number and, in fact, observed that neglecting q 2 slightly 
improved the agreement between the LES and DNS 
results. 
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Vreman et al 9 confirmed the above results with 
their simulations of a three-dimensional temporal 
compressible mixing layer at a mean convective Mach 
number of^0.2. In a priori tests, the SGS model that 
neglects q *" was found to give a better correlation with 
the DNS results. Furthermore, LESes conducted with 
the dynamic SGS model for q~ were unstable for the 
cases that were studied. 

For low-Mach number turbulence LES, neglecting 
the term q will not introduce large errors in the results 
and is actually desirable in some cases, as the above 
findings showed. As a result, the term q is neglected 
in this study. This assumption is analogous to the 
Stokes assumption for the viscous stress^tensor in the 
Navier-Stokes equations. With the term q' omitted, the 
SGS stresses can be included in the Navier-Stokes 
equations by simply replacing the laminar viscosity 
coefficient )i with \i € ^ where = P + Psgs anc * 
Uses = CpA 2 |s| . 

Lar ge-Eddv Simulation of 
Turbulent Flow in a Square Duct 

To validate the numerical method for turbulence 
simulations of duct flows, LES of fully developed 
turbulent flow in a square duct is performed. For the 
purpose of comparison, a low-Reynolds number, square 
duct DNS solution is available. This DNS database was 
used by Mompean et al. 19 to evaluate nonlinear k-e 
turbulence models. Another DNS solution of the fully 
developed turbulent square duct flow at a slightly lower 
Reynolds number is also available. 20 

Figure 4 shows the coordinate system and geometry 
for the square duct flow. In this test case, the Reynolds 
number based on the mean streamwise velocity and 
hydraulic diameter is 4800. Based on the friction 
velocity and hydraulic diameter, the Reynolds number 
is 320. 

Table 2 shows a summary of the flow properties for 

the test case, assuming an average Mach number of 0.3 

and standard sea level properties for air. The 

computational domain size used in the LES is 

\2D h x D h x D h . In choosing the size of the 

computational domain, care must be taken to ensure that 

the length of the computation domain is large enough to 

adequately contain the largest turbulence structure. 

Two-point velocity correlations for three different 

cross-stream positions in the duct were computed from 
r on 

the DNS solution by Gavrilakis The correlations for 
all three velocity components become essentially 0 at a 



Figure 4. Coordinate system and geometry for the 
square duct flow. 

duct length of approximately 6 D H , so that a length of 
12 D h should be adequate to capture the streamwise 
turbulence structures. Two different grids are used in the 
present LES, and table 3 shows the simulation 
parameters. The sampling time for the turbulent 
statistics is large compared to the time step size, but is 
small compared to the eddy turnover time in order to 
capture the unsteadiness of the turbulent flow. 


Table 2. Flow properties for the turbulent square duct 
test case. 


F low properties 

Values 

Average Mach number, 

0.3 

Average st'eamwise velocity, u ave 

102.4 m/sec 

Average friction velocity, u x 

6.83 m/sec 

D P u ave D H 

R e = 

p 

4800 

P M t d h 
We. = 

1 H 

320 

Hydr tulic diameter, D H 

7.37 x 10 -4 m 

Mean \ ressure gradient, P 

-289,646 Pa/m 

0.5 D h 

Eddy turnover time, 

“t 

5.4 x 10- 5 sec 
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Table 3. Parameters for the LES 


Parameters 

Grid A 

Grid B 

Grid size (I x J x K) 

129x90x90 

257 x 100 x 100 

Number of cells 

1,013,888 

2,509,056 

Minimum 
resolution 
(in wall units) 

30 x 1.88 x 1.88 

15 x 1.69 x 1.69 

Maximum 
resolution 
(in wall units) 

30x4.86x4.86 

15x4.37x4.37 

Sampling time, 
A t s (sec) 

6.0 x 10“ 7 

5.0 x 10’ 7 

Time step size, 
At (sec) 

3.0 x 10~ 9 

2.5 x 10~ 9 

CFL number 

0.98 

0.93 


The periodic boundary condition used for the inflow 
and exit boundary of the square duct is similar to the one 
used by Coleman et al. 21 With this boundary condition, 
all of the flow conditions are periodic at the duct inlet 
and exit planes. The driving pressure gradient in the 
duct is specified in the flow equations as an extra body 
force term. 

To reduce the number of iterations required for 
convergence, the initial conditions for the large eddy 
simulations were obtained from interpolating a DNS 
solution provided by Gavrilakis. The DNS was done 
using a 768 x 127 x 127 grid. The current simulations 
were done using two different grid sizes, 129 x 90 x 90 
(grid A) and 257 x 100 x 100 (grid B). The finer LES 
grid B, which gave good results, is approximately 
20 percent of the total size of the DNS grid. 

The convergence of the LES is determined by 
monitoring the time history of the total wall shear stress. 
For fully developed turbulent flow in a straight square 
duct, conservation of the mean streamwise momentum 
shows that the mean driving pressure gradient and the 
total wall shear stress are related by the following: 

j X w dA = -VP g (49) 

The surface integral is over the four side walls of the 
square duct, so that A s = 4x12 x D H 2 . P g is the 


mean driving pressure gradient, and V is the total 
volume of the duct, given by \2x D H *. Defining 

X, = — fx.cM and P n = the familiar relation 
H A* >v * dx 

between the mean pressure gradient and wall shear 
stress in fully developed flow in a square duct can be 
recovered. 



Figure 5 shows the time history of the total side wall 
shear force for the LES using grid B. The instantaneous 
side wall shear force level from the LES, shown as a 
solid line, is seen to fluctuate about a mean value, 
indicating that flow equilibrium has been reached in the 
current simulation. The time average of the computed 
side wall shear force is 0.001387 N. This value is in 
excellent agreement with the exact value of 0.001 391 N, 
computed from equation (50) and shown as the 
horizontal dashed line (fig. 5). Because the time step is 
constant, the number of iterations shown in figure 5 is 
directly proportional to the time elapsed. The simulation 
was conducted for 218,600 time iterations, which is 
approximately 10 eddy turnover times (as defined in 
table 2). 



Number of iterations (in thousands) 

980503 

Figure 5. Convergence history for the total wall shear 
stress. 

The parallel implementation and the results of the 
parallel performance studies have previously been 
published. 1, 3 The code was implemented on parallel 
computer systems using the message-passing 
programming model and message-passing libraries such 
as Message-Passing Interface (MPI) and Parallel Virtual 
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Machine (PVM). The parallel speedup was found to be 
very good, especially for large numbers of grid points/ 1 
Using 128 processors on the T3D computer (Cray 
Research, Inc.; Eagan, Minnesota), the simulation of 
these 10 eddy turnover times (5.5 x 10 -4 sec in physical 
flow time) took 772 hr or approximately 1 month of 
central processing unit time. The same simulation 
would have taken approximately 150 hr on an SP2 (IBM 
Corporation, Austin, Texas) or T3E (Cray Research, 
Inc.) computer. Regardless of the computer platform 
used, this computational time is a large cost and shows 
that even with the parallel computer systems available 
today, turbulence simulation is still a formidable task. 
However, parallel CFD algorithms that can efficiently 
scale up with extremely large numbers of processors 
offer the only real hope that turbulence simulations can 
be done in a reasonable amount of time in the future. 

The LES results shown in figures 5 to 14 are obtained 
using grid B and a modified Roe FDS with an £j value 
of 0.03. The modification to the Roe FDS and the 
definition of the £j parameter will be discussed below. 

Figure 6 shows the mean streamwise velocity profile 
along a wall bisector. The LES solution (solid line) is 
compared with the DNS solution (diamond dots) 
supplied by Gavrilakis. The mean velocity profile in the 
LES was averaged both in time and space. The 
agreement can be seen to be very good. 


LES 



U ave 980504 


Figure 6. Mean streamwise velocity profile for fully 
developed turbulent flow in a square duct along the wall 
bisector. 

Figure 7 shows the mean secondary velocity vectors 
from the LES. In straight ducts of noncircular cross- 
sections, turbulence-driven secondary flows are known 
to exist. These flows are different from the pressure- 
driven secondary flows found in curved ducts. In 



0 .5 1.0 

*/°H 980505 


Figure 7. Mean secondary velocity vectors from LES 
with a 257 y 100 x 100 grid. 


straight square ducts, the turbulence-driven secondary 
flows are directed from the center of the duct toward the 
comers along the comer bisectors, and have been found 
to be produced by the anisotropy of the Reynolds 
stresses in the cross-sectional plane of the square duct. 22 
Although tfe magnitudes of these secondary velocities 
are extremely small compared to the mean average 
streamwise velocity (on the order of 2 percent in this 
simulation) these velocities have been found to be 
important features of this flow. 

Figure 7 shows that the comer vortices produced by 
the secondary flows are captured in this simulation. 
Although some asymmetry is still evident in the plot, 
the overall features of the secondary flows are well- 
predicted b; the current simulation. 

To deter nine the accuracy of the simulation in 
capturing turbulence-driven secondary flows, the mean 
secondary velocity profiles along the lines 
z/(0.5 D H ) = 0.15, 0.30, 0.50, 0.70, and 0.80 are 
compared with the DNS results in figures 8 to 12. 
Generally, $ ood agreement is obtained between the LES 
and DNS m^an secondary velocity profiles. 

Figures 13 and 14 show the turbulence statistics. In 
figure 13, tie mean Reynolds stress profile along a wall 
bisector is compared with the DNS solution, and in 
figure 14, tie mean turbulence intensities w rms , v rms , 
and vv rms a e plotted. These results have been quadrant- 
averaged as well as averaged in space and time. 
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LES 

o DNS 




Figure 8. Mean secondary velocity profiles along 
z/(0.5 D h ) = 0.15. 


Although the agreements between the LES and DNS 
solutions are seen to be very good for this simulation, 
the accuracy of the LES in capturing the turbulence 
velocity fluctuations was found during the turbulence 
simulations to be highly dependent on the numerical 
dissipation and the grid size used. The effects of the Roe 
FDS upwinding term and grid size on the computed 
turbulence velocity fluctuations are examined next. 

Effect of the Roe 
Flux-Difference Splitting Term 

Although the Roe FDS implemented in this code gave 
good results for Euler 2 and laminar Navier-Stokes 3 test 
cases, the full Roe FDS term was found to be too 
dissipative for LES. Incorrect levels of turbulent 
velocity fluctuations are obtained when the normal Roe 
FDS term is used in turbulence simulations. This 
problem was solved with a simple modification to the 
Roe FDS algorithm. In equation (10), the inviscid fluxes 
normal to a cell boundary is approximated as 


LES 

o DNS 



W / U T 980507 


Figure 9. Mean secondary velocity profiles along 
z/(0.5 D h ) = 0.3). 

f = ^l + ^r)-^I(Ur-U l ) 

This approximation can be interpreted to state that the 
normal component of the inviscid flux at a cell boundary 
is the sum of the central difference of the fluxes on the 

left and right states, ^(f L + f R ), P^ us R° e 

upwinding dissipation term, -j-|A|(U R - U L ) j . If this 

interpretation is used, then the amount of Roe 
upwinding dissipation can be controlled using a 
multiplying factor in front of the Roe FDS term, such 
that 

f = ^ f L + f R)-e.{^(U R -U L )} (51) 

where £j can range between 0 and 1 . e | = 0 corresponds 
to central differencing only, and £j = 1 corresponds to 
the full Roe FDS. 
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Figure 12. Mean secondary velocity profiles along 
z/( 0.5 D h ) = 0.8). 


LES 

* DNS 

1.0 
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_ .6 

- u V/u* 
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0 .2 .4 .6 .8 1.0 

y/0‘ 5D H 980511 

Figure 13. Mean Reynolds stress profile along the wall 
normal bisector. 





Y/0*5D H 980512 


Figure 14. Turbulent intensities along the wall normal 
bisector. 


Note that Lin et al., using the same interpretation of 
the Roe upwinding term as equation (51), also 
concluded that the normal Roe upwinding term 
produces too much numerical dissipation for 

computational aeroacoustics applications. 2 ' Lin et al. 
found that using e, values of approximately 0.1 gave 
good results for acoustics computations."' 

In the LES conducted here, values of less than 0.1 
are needed to give good turbulence results. Omitting the 
Roe FDS term altogether (£j = 0) causes all calculations 
to be unstable, and the best turbulence solutions are 
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obtained using the smallest possible values of e, that 
can still provide stable calculations. In general, the finer 
the grids, the smaller the minimum values of E ] that can 
be used. For a grid of 129 x 90 x 90, the minimum value 
of for stable calculation is 0.05, and for the 
257 x 100 x 100 grid, the minimum value is 0.03. The 
LES results presented in the preceding section were 
obtained using £j = 0.03. 

To study the effect of the Roe upwinding term on the 
turbulence simulation, LES are made for the same grid 
size of 129 x 90 x 90 but with different values of £j. 
Figure 15 shows the effect of Roe FDS on the mean 
streamwise velocity profile. Near the wall, using the full 
Roe FDS term produces a mean velocity gradient that is 
much less than both the DNS solution and the LES 
solution with the reduced Roe FDS. A similar effect is 
observed in figure 16, where the mean Reynolds stress 
profile obtained with £j = 1.0 is much lower than 
expected. 

Figure 17 also shows the excessive numerical 
dissipation of Roe FDS in the turbulence solution. In 
this figure, the solution with £j = 1.0 gives a 
significantly higher level of w rms and lower levels of 
v fms and tv . Although they did not use the Roe FDS, 
Wang and Pletcher 24 reported the same problem in their 
LES of fully developed turbulent channel flow using an 
upwind CFD algorithm. 

From studying the results shown in figures 15 to 17, 
the full Roe FDS upwinding dissipation can be seen to 
be detrimental to the turbulence solution, and an 


LES with £<| =1.00 

LES with £ 1 = 0.05 

o DNS 



Figure 15. Effect of Roe FDS on the mean streamwise 
velocity profile. 


LES with £ 1 = 1 .00 

LES with = 0.05 

o DNS 



Figure 16. Effect of Roe FDS on the mean Reynolds 
stress profile. 


improvement in the solution quality can be obtained by 
reducing the contribution of the Roe upwinding term. 
However, continuing to reduce the contribution of the 
Roe term until a good agreement is achieved is not 
possible. For a given grid size, a minimum amount of 
Roe FDS upwinding dissipation is required for stability. 
For the 129 x 90 x 90 grid, the minimum £ ( value for 
numerical stability is 0.05, and values smaller than this 
minimum will cause the calculation to be unstable. To 
improve the accuracy of the turbulence simulation, 
using a finer grid that in turn allows a smaller £ } value 
to be used is necessary. In the next section, the effect of 
a finer grid on the quality of the turbulence solution will 
be studied. 

Effect of the Grid Size 

The previous section showed that reducing the 
contribution of the Roe FDS term will improve the 
quality of t le solution. But using a 129 x 90 x 90 grid, 
reducing £ to the minimum value of 0.05 still does not 
give a good agreement with the DNS solution. In this 
section, a filer grid with 257 x 100 x 100 points is used, 
and the minimum value of £j that can be used is 
lowered to ).03. 

Figure 8 shows a comparison of the mean 
streamwise velocity profiles. Using the finer grid in the 
LES produces an almost perfect agreement with the 
DNS solution. Figures 19 and 20 show the same 
improveme it in the profiles of the turbulence intensities 
and the mean Reynolds stress, respectively. 
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LES with £ 1 = 0.05 

o DNS 
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Figure 17. Effect of Roe FDS on the turbulent 
intensities. 


Another effect of grid size can be observed by 
examining the mean secondary velocity vectors in a 
cross section of the square duct. Figure 21 shows the 
solution obtained using grid A and e 1 = 0.05. The 
secondary turbulence-induced velocity field is captured 
by the coarser grid. However, comparing this result with 
the finer grid B result in figure 7 shows that the corner 
vortices in figure 7 are somewhat smaller than those in 
figure 21. Figure 22 shows the streamwise-averaged 


LES with 129 x 90 x 90 grid 

LES with 257x 1 00x1 00 grid 

o DNS 



Figure 18. Effect of grid size on the mean streamwise 
velocity profile. 


instantaneous secondary velocity vectors from the DNS 
solution for reference. For the DNS grid, the near wall 
vortical structures are even smaller than either of the 
LES grids. This comparison shows that the near wall 
turbulent structures are better resolved with finer grids. 

Effect of the Subgrid-Scale Model 

The previous section showed that the LES solution 
with the fine grid gives the best agreement with the DNS 
solution. Although the LES fine grid is only 20 percent 
of the size of the DNS grid, the LES grid density in the 
crossflow plane is 60 percent of the DNS crossflow 
plane density. As the LES grid resolution in the 
crossflow plane approaches the DNS resolution, the 
effect of the SGS model on the turbulence solution for 
this particular LES grid is interesting to see. An 
additional simulation was performed using grid B, the 
finer LES grid, with no SGS model. This simulation is 
effectively a coarse grid DNS. Figure 23 shows a 
comparison of the mean streamwise velocity profiles. 
The use of the SGS model makes essentially no 
difference in the mean streamwise velocity solution. 
Small differences are also observed in the turbulent 
velocity fluctuations shown in figure 24. Figure 25 
shows a comparison of the mean Reynolds stress 
profiles. The biggest difference is at the peak of the 
Reynolds stress profile, where the LES solution with no 
SGS model predicts a slightly higher peak than the 
DNS, and the LES solution with the SGS model predicts 
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Figure 22. Stream wise-averaged secondary velocity 
vectors from the DNS solution. 
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Figure 23. Effect of the SGS model on the mean Figure 24. Effect of the SGS model on the turbulent 
streamwise velocity profile. intensities. 


a slightly lower peak than the DNS. Note that the LES 
turbulence statistics were only computed for the 
resolved velocity field. As a result, the LES Reynolds 
stress profile with the SGS model should be lower than 
the DNS solution because the SGS contribution was 
not included. 


These results indicate that an SGS model is not 
needed for an accurate simulation of this test case. As 
discussed earlier, the grid resolution in the near wall 
region has to be very fine to resolve the small energy- 
producing structures there. The fine grid LES conducted 
here has effectively approached the DNS in the near 
wall limit. 
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Figure 25. Effect of the SGS model on the mean 
Reynolds stress profile. 


Conclusion 

A new, parallel, finite-volume computational fluid 
dynamics algorithm was developed for large-eddy 
simulation (LES) of turbulent flows using parallel 
computer systems. Major components of the algorithm 
included piecewise linear least-square reconstruction of 
the unknown variables, trilinear finite-element 
interpolation for the spatial coordinates, Roe flux- 
difference splitting (FDS), and second-order 
MacCormack explicit time marching. The parallel 
implementation was accomplished using the message- 
passing programming model. 

For the first time, a parallel, unstructured, finite- 
volume numerical algorithm was used for LES of 
turbulent flow in a square duct, and several conclusions 
have been drawn regarding the accuracy and efficiency 
of this numerical algorithm. Comparison with the direct 
numerical simulation (DNS) solution showed that the 
standard Roe FDS upwind dissipation adversely affects 
the accuracy of the turbulence simulations. A 
modification to the standard Roe FDS method was 
proposed in which the inviscid flux is computed as the 
arithmetic average of the right and left fluxes plus 
the product of the Roe FDS dissipation term and a 
reduction factor. For accurate turbulence simulations, 
only 3-5 percent of the normal Roe FDS dissipation 
was found to be needed. 

The finer, 257 x 100 x 100 LES grid required less Roe 
FDS upwind dissipation for stability and produced a 
more accurate solution than the 129 x 90 x 90 LES grid. 
The near wall vortical structures were better simulated 


by the finer grid LES, and the effect of the subgrid-scale 
model on tie accuracy of the results was found to be 
small for th^ fine grid LES, which is nearly as fine as the 
DNS grid in the near wall region. 
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